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/~s _ Abstract. I review recent numerical studies of accretion disks, focusing on measure- 

■ ments of the turbulent shear stress, or a, in the shearing box model. I conclude with a 
I list of astronomically relevant open questions that can be settled via future numerical 

experiments. 
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S : INTRODUCTION 
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The last decade has brought a flood of cheap, fast workstations and easy access to 

■ supercomputers to the astronomical community. For those struggling to understand 
complicated phenomena like the active nucleus of NGC 4258, with its beautifully 
precise maser spots orbiting in a thin, warped disk, these machines seem to offer 
the seductive possibility of creating a detailed and comprehensive simulation. One 
could imagine that by including all the relevant physics, such a simulation would 

', self-consistently generate X-rays from the relativistic accretion flow in the neigh- 

^ ' borhood of the black hole (ending the debate between proponents of ADAF and 

disk corona models); it would warp straightaway (giving away the origin of the 
^ . warp); it would include a chemical network to calculate the local abundance and 

excitation of water; it would perform three dimensional (3D) radiative transfer in 
the water maser lines; and of course, it would produce a jet. Our curiosity would 
be completely satisfied. 

Of course, such global and physically complete simulations are mere fantasy. 
They are not possible now, nor will they be except in some distant, post-Moore's- 
law future world. It is not even clear that they would be desirable! Even if one could 
determine the right initial conditions — probably impossible even in principle- the 
output would be so complicated (cf. [1] on spiral structure) that it would not be 
clear what was important. 

Absent global, ab initio simulations one is reduced to solving simplified model 
equations that rely on reasonable but untested physical assumptions. Thus models 
of dwarf novae (see, e.g. [2,3]) typically evolve an azimuthally averaged and height- 
integrated set of equations for the surface density E: 



O 
H 



dtE{r, t) = \ {-^dr{r^^r^) - ^) - tw. (1) 

Here S = surface density, Vl = rotation frequency, I^w = mass lost per unit area in 
winds, Wr^ = J dz Wrcf, = height-integrated shear stress, and r = direct torque per 
unit area- possibly supplied by a magnetohydrodynamic (MHD) wind or by the 
tidal field of the secondary. Equation (1) is in a sense fundamental; it is derived 
from the angular momentum and continuity equations in the limit that the disk is 
thin. 

In solving equation (1) two assumptions are almost always made. First, r = 
Evt^ = 0. This assumption is required by our ignorance of the very difficult, global 
problem of disk winds; the relative importance of external torques r and internal 
stresses Wrtj, remains one of the outstanding problems of disk physics. Second, 
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Wr<l, = -P^y; y = OL(?J^, (2) 

where = sound speed, v is the "anomalous viscosity" and a its nondimensional 
counterpart. This equation is not fundamental; it is the simplest possible represen- 
tation of the effects of turbulence with the correct dimensional form. 

So perhaps if we cannot, god-like, simulate an entire active galactic nucleus on the 
computer, we can at least set the more modest goal of using numerical experiments 
to understand the origin and evolution of w^^. 

The turbulent shear stress, or "anomalous viscosity," has accreted an aura of 
mystery over the years that is not entirely deserved. We know what the basic 
governing equations are. They do not involve exotic particle physics. In many 
cases (ADAFs and disk coronae are potentially important exceptions) they do not 
even involve plasma kinetics. One may simply treat the disk as a magnetized fluid. 
Then 
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where = gravitational potential. Any change in the angular momentum of the 
fluid in the disk is due to a torque N = r x Dw /Dt. Torques can only be due, then, 
to pressure gradients, magnetic forces, or nonaxisymmetric gravitational fields. Ra- 
diation forces, which have been dropped here, are negligible except in certain special 
circumstances. 

The fluid equations are easy to write down but hard to solve. To make a numerical 
solution practical, one wants to (1) include the minimal relevant physics, (2) take 
advantage of a symmetry, such as axisymmetry, if possible, and (3) start with a 
simple disk model. 

The first of these considerations motivated early workers on disk dynamics to 

neglect magnetic forces in the equation of motion. Since the discovery by Balbus 
& Hawley [4] of a linear instability in weakly magnetized disks, it has become clear 



that this approximation, while nobly motivated, misses perhaps the most important 
physics (some disks, however, may be so neutral that they are decoupled from the 
magnetic field [5,6]). 

One might still hope to use axisymmetry. Alas, early axisymmetric models of 
magnetized disks showed the development of structures called "channel solutions," 
consisting of superposed layers of material with large radial velocities [7]. An 
analytical study [8] showed that these structures are unstable in three dimensions. 
The upshot is that two dimensional experiments have a nonlinear outcome that 
is completely different from that of three dimensional experiments. While it is 
possible that channel solutions may be relevant under certain special conditions, 
three dimensional experiments are really required to advance our understanding of 
disks. 

Finally, one wants to start with a simple disk model. A natural choice is the 
"local model," which is a rigorous first order expansion of the equations of motion 
in e = H{r)/r <^ 1 {H{r) = disk scale height) in a frame comoving with the disk. 
Boundary conditions are also needed so that the local model can be mapped to a 
finite computational domain. The "shearing box" boundary conditions [9-11] are 
well suited to this. They are similar to periodic boundary conditions except that 
they allow for the presence of shear due to differential rotation. They permit the 
study of a small piece of the disk without reference to poorly understood inner 
and outer radial boundaries, although vertical boundary conditions must still be 
supplied. 

All this motivates interest in numerical experiments in the shearing box. 1 will 
review recent work, then discuss open questions that might be addressed with more 
work, and more CPU cycles, in the future. 



NUMERICAL EXPERIMENTS AND RESULTS 

The 3D shearing box experiments done to date have all considered a compressible 
fiuid using finite difference or finite-difference-like methods ^ They have included a 
variety of physics: magnetic fields and pure hydrodynamics; resistive [12] and am- 
bipolar [14] diffusion; other-than-Keplerian rotation curves [15]; forced convection 
from heating at the midplanc [16,17]; and self-gravity [18]. 

Shearing box experiments may be divided into two broad classes. The first class 
are unstratified, i.e. neglect the vertical structure of the disk. Thus if is the 
gravitational potential, d(j)/dz — and so in the initial laminar equilibrium dp/dz — 
0. The second class are stratified. They include the usual vertical tidal potential 
(f) — \^'^z^, so an isothermal equilibrium would have p — poexp{—z'^/{2H'^)). 



A incompressible study of MHD turbulence in disks using pseudo-spectral methods is also 
possible and would provide an interesting check on the numerics. 



Unstratified, Magnetized Shearing Box 



Why not dispense with unstratified models and go straight to the more reahstic, 
stratified case? First, the unstratified models are numerically less demanding. The 
stratified boxes include low density regions that require small timesteps because 
of the Courant condition in a magnetized fluid. Stratified boxes also assign only 
a fraction of the zones to the turbulent center of the disk; the rest are assigned 
to the nearly force-free disk atmosphere. Second, the unstratified boxes allow the 
study of the nonlinear development of the Balbus-Hawley instability independent 
of buoyancy effects. A final, post hoc justification is that stratification is not yet 
measured to have a significant effect on the outcome. This may, however, not be 
true in future large, highly resolved experiments. 

The model has several important dimensionless parameters. One set of param- 
eters describes the shape and size of the box relative to H — c^/Q. Typically 
Lr X Lff) X Lz = {1 X 271 X 1)H. It is natural to set = (1 — 2)H; this simulates 
the limited vertical scale available in a disk. The other dimensions are limited only 
by CPU time. Another set of dimensionless parameters describes the numerical 
resolution. Most experiments to date use 2^ to 2^ zones along each axis. A final 
set of parameters describes the mean magnetic field (B) (the brackets denote a 
spatial average), made dimensionless by comparison of the Alfven speed with Cg. 
The boundary conditions force (Bz) and {Bj.) to be constant in time, while (B^) is 
fixed if and only if {Br) = 0. 

Two groups have reported unstratified shearing box experiments [11,19,12]. All 
show the development of MHD turbulence initiated by the Balbus-Hawley insta- 
bility; the turbulence eventually settles down into a final state that is independent 
of all the initial conditions except the mean magnetic field. 

There are many interesting quantities that can be measured in the final state. I 
will focus on just one: the shear stress in the nondimensional form 



where I have assumed the disk is Keplerian. The first term is the magnetic, or 
"Maxwell," stress, and the second term is the fluid, or "Reynolds" stress. The 
brackets indicate an average in r,(f),z and t. 

The single most important result to emerge from the unstratified shearing box 
experiments is that the Balbus-Hawley instability initiates MHD turbulence that 
has a ^ 10~^. No other proposed mechanism for angular momentum transport 
(except gravitational instability) has produced this result in a purely local disk 
model. 

Measured values of a depend on the box size, on the mean magnetic field strength, 
and at worst very weakly on numerical resolution [11,12]. Crudely speaking. 
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where (V^) = mean Alfven velocity. Thus for weak mean fields, a ~ 10~^, while 
intermediate strength mean fields raise a above this base level. For mean fields 

that are strong in the sense that j is large compared to LiVt the field is linearly 
stable and equation (5) docs not apply (sec [4,13] for linear stability criteria). The 
flow is laminar, and a — 0. Equation (5) should not be taken too seriously as yet; 
while the sense is likely correct, the coefficients could change substantially in future 
experiments. 

Notice that in the absence of a mean field all the currents that sustain the field 
are contained within the computational domain and so are subject to decay. In this 
sense, the zero-mean-field experiments demonstrate the existence of a dynamo: they 
show that MHD turbulence in disks, driven by the Balbus-Hawley instabihty, can 
sustain a magnetic field in the presence of dissipation [12]. 

Finally, some zero-mean-field experiments [12] included a finite resistive diffu- 
sivity Tj. These show that a is sensitive to the presence of resistive diffusion. 
Defining the magnetic Reynolds number Rcm = CgH/r), turbulence decays at 
ReM ~ 10^ — 10^. Dwarf novae disks in quiescence may fall in or below this 
range [6]. 

Magnetized, Stratified Shearing Box 

Stratified shearing box experiments have been carried out by two groups 
[14,20,21,15]. The experiments have = 4-6 X\/2H, and numerical resolution 
similar to the unstratified models. Both groups report generally consistent results, 
which is remarkable since very different numerical methods were used. The strati- 
fied experiments give a ~ 10^^'^. 

The stratified boxes produce a slightly lower than the unstratified boxes. Part 
or all of this difference may be due to the lower effective numerical resolution of 
the more numerically demanding stratified experiments. The stratified boxes are 
sensitive to resolution in that, when the resolution is increased, a also increases 
[21,20]. Once higher numerical resolution becomes practical, the stratified experi- 
ments ought to be repeated until convergence can be demonstrated. 

The measured a is insensitive to the vertical boundary conditions. Again, this 
may change in larger, converged experiments. Some experiments have used vertical 
boundary conditions that allow the mean field to evolve [14,20,15]. In my view this 
is risky because the mean field is generated by currents "in the boundary" and not 
within the computational volume itself; changes in the mean field are thus generated 
by interaction of the fluid with the boundary conditions. This is, however, mainly 
a matter of taste. It is fair to say that we do not yet know what the most relevant 
vertical boundary conditions are. 

One of the main motivations for the stratified experiments is that they can in 
principle be used to estimate the vertical run of turbulent dissipation in disks, 
thereby removing a serious obstacle to predictive models for disk spectra. In isother- 
mal models the run of dissipation is not calculated directly, but two closely related 



quantities are: Wr(j,{z) and Sz{z), the vertical component of the Poynting flux. It 
is found that [21] 

S,{z = //) ~ 0.01 dz InwrJz). (6) 
Jq 2 

Thus only a small fraction of the power extracted by the turbulent shear stress from 
the differential rotation emerges as MHD waves. Larger values of a in future exper- 
iments, however, would imply stronger fields and hence greater magnetic buoyancy. 
It is also found [21] that a is not constant with height, as is commonly assumed. 
A better fit to the data is given by a ~ (p/po)~^^^, consistent with equation (5). 
Clearly the situation with the stratified experiments is not entirely satisfactory; it 
is a challenge for future experiments to make the dissipation and vertical energy 
transport explicit in a converged calculation. 

Unmagnetized Box 

Are there any local transport processes that can compete with or dominate MHD 
turbulence? And what transport processes govern the evolution of disks that are 
nearly neutral and thus poorly coupled to the magnetic field? 

Convectively driven turbulence was once thought a promising transport mecha- 
nism [22,23]. There was an early warning that something might be amiss, however, 
from a quasilinear study [24] which showed that nonaxisymmetric convective modes 
produce an inward angular momentum flux. Subsequent shearing box experiments 
[16,17] showed that both forced convection and overturning of an initially unstably 
stratified disk led to inward angular momentum transport [a < 0) in the fully 
nonlinear regime. This counterintuitive result is a nice illustration of the value of 
numerical experiments. 

There is still one regime in which convectively driven angular momentum trans- 
port could play a role. That is in geometrically thick fiows with unstable radial 
stratification- for example ADAFs (see [25] for a review). A shearing box experi- 
ment with forced radial convection might be revealing, but the shearing box is not 
a good model for these flows. Global models are really required. 

Shearing box experiments have also permitted the direct evaluation of another 
once-promising mechanism for initiating turbulence in disks: nonlinear hydrody- 
namic instability [26]. In these experiments a purely hydrodynamical Keplerian 
shear flow is violently perturbed. It is found that the flow returns to a laminar 
state, independent of the amplitude of the initial perturbation. 

Because the numerical experiments are run at a Reynolds number that is low 
in comparison to that in astrophysical disks, nonlinear hydrodynamic instability 
cannot be rigorously ruled out. But analytic arguments [26] (see also Balbus, this 
volume) and the absence of nonlinear instability in laboratory analogs relevant to 
Keplerian shear flow (see [27] and references therein) make it likely, in my view, 
that Keplerian disks are nonlinearly hydrodynamically stable. 



Finally, self-gravity may be noticeable or even dominant in disks around young 
stars and in active galactic nuclei. The shearing box has been used to study grav- 
itational instability in cooling disks [18]. Cooling is essential since in its absence 
the disk simply heats up until it is stable. The outcome is a fluctuating state with 
Toomre's stability parameter {Q) ~ 1 and significant outward transport of angular 
momentum via gravitational and Reynolds stresses. These experiments also show 
that self-gravity produces truly local transport in a sense to be made clear below. 

OPEN QUESTIONS 

I will conclude with a short list of open questions that are astronomically relevant 
and can be answered in the near future with numerical experiments in the shearing 
box. 

1. Is a determined locally? A different way of phrasing this question is, does a 
converge as the planar box size L,. and arc increased? Experiments to date find 
turbulence with most of the energy in structures with scale comparable to the box 
size. Thus the outcome is limited by the box size. 

If a is determined locally one expects that the autocorrelation function of fluid 
variables to decay rapidly on scales > H. Equivalcntly, fluid variable power spectra 
should turn over and decline at small fc,. and k^. Because the natural scale for 
the turnover is H, it should be seen in boxes only slightly larger than current 
experiments. The detection of the turnover would be a significant milestone for 
disk physics in that it would show that a is locally determined and so validate the 
use of the local model. 

The alternative is that the fiuid autocorrelation functions decline only slowly or 
not at all with distance in the plane of the disk. Then the shear stress at any point 
in the disk can be influenced by conditions at points many scale heights away. This 
is not inconceivable. The a ~ [H/rY prescription common in studies of dwarf 
novae, for example, requires this sort of nonlocality because local turbulence must 
somehow "know" about the large scale structure of the disk and in particular the 
local radius r. If a is nonlocal, however, we must abandon the shearing box and 
move to global models. In what follows I assume that a is local and thus that the 
shearing box is relevant. 

2. How does a depend on resistive diffusion rj, viscosity v, and their ratio, the 
magnetic Prandtl number Pru = I'/rj? In most circumstances both resistivity and 
viscosity are negligible, but even then, it has been argued, Pru may govern the 
character of MHD turbulence (see [27] and references therein). This hypothesis 
can be tested with sufficiently high resolution experiments that allow a reasonable 
separation of the resistive and viscous lengthscales. Resistivity is not always neg- 
ligible, however. Recent work [6] shows that the standard disk instability model 
for dwarf novae implies a magnetic Reynolds number of order 10^ for dwarf nova 
disks in quiescence. Such low Rcm may have a direct impact on the development of 
MHD turbulence (cf. [12]). Protostellar disks also suffer from high resistivity; some 



parts are likely to be completely decoupled from the magnetic field [5] . This issue 
has direct astronomical relevance and can be studied with codes and computers 
that are now available. 

3. How does a vary in time? Local model MHD experiments have so far only 
sought to measure a when the disk scale height is steady in time. Some of the most 
potentially revealing phenomena in disk systems, however, involve rapid changes in 
disk temperature. Under these circumstances, how long does it take a to readjust 
to its steady state value? There is already preliminary evidence that this relax- 
ation time is long in that the time required for the zero-mean-field unstratified 
experiments to reach a steady state is many orbital periods [12]. But a more direct 
determination can be made in stratified, isothermal shearing box experiments in 
which the temperature is forced to change suddenly. If the relaxation time is long 
compared to the thermal timescale, then the structure of cooling and heating fronts 
in dwarf nova disks may be quite different from what is now imagined. 

4. How does a vary near disk edges? Gaps may be opened in disks around 
young stellar objects by stellar or planetary companions. The size of the gap and 
accompanying torque on the companion are sensitive to the surface density profile 
close to the disk edge, which is determined by a balance between tidal torques on 
the disk and all other torques. One way of studying disk edges would be to allow 
a point mass to orbit in the shearing box and clear a gap in the disk. Such an 
experiment could test whether a varies near the edge of the disk in a way that 
is consistent with the usual treatment of Wnj, as a viscous stress. Perhaps it does 
not; perhaps a increases rapidly within a few x2'nH of the edge. Such experiments 
could also test whether magnetic fields drive cross-gap accretion. Numerical work 
to date has focused on the nonmagnetic problem (e.g. [28]), but this approximation 
may be missing an important piece of the puzzle. 
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